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We introduce two- and one-dimensional (2D and ID) models of a binary BEC (Bose-Einstein 
condensate) in a periodic potential, with repulsive interactions. We chiefly consider the most funda- 
mental case of the inter-species repulsion with zero intra-species interactions. The same system may 
also model a mixture of two mutually repulsive fermionic species. Existence and stability regions 
for gap sohtons (GSs) supported by the interplay of the inter-species repulsion and periodic poten- 
tial are identified. Two-component GSs are constructed by means of the variational approximation 
(VA) and in a numerical form. The VA provides accurate description for the GS which is a bound 
state of two tightly-bound components, each essentially trapped in one cell of the periodic potential. 
GSs of this type dominate in the case of mtra-gap solitons, with both components belonging to 
the first finite bandgap of the linear spectrum (only this type of solitons is possible in a weak lat- 
tice). Inter-gap solitons, with one component residing in the second bandgap, and intra-gap solitons 
which have both components in the second gap, are possible in a deeper periodic potential, with 
the strength essentially exceeding the recoil energy of the atoms. Inter-gap solitons are, typically, 
bound states of one tightly- and one loosely-bound components. In this case, results are obtained 
in a numerical form. The number of atoms in experimentally relevant situations is estimated to 
be ~ 5,000 in 2D intra-gap soliton, and ~ 25,000 in its inter-gap counterpart; in ID solitons, it 
may be up to IC. For 2D solitons, the stability is identified in direct simulations, while in the ID 
case it is done via eigenfrequencies of small perturbations, and then verified by simulations. In the 
latter case, if the intra-gap soliton in the first bandgap is weakly unstable, it evolves into a stable 
breather, while unstable solitons of other types (in particular, inter-gap solitons) get completely 
destroyed. The intra-gap 2D solitons in the first bandgap are less robust, and in some cases they are 
completely destroyed by the instability. Addition of intra-species repulsion to the repulsion between 
the components leads to further stabilization of the GSs. 

PACS numbers: 03.75.Lm, 05.45.Yv 



I. INTRODUCTION 



Solitons in Bose-Einstein condensates (BECs) have drawn a great deal of attention as robust nonlin- 
ear matter-wave pulses. Bright (localized) solitons were first experimentally created in an effectively one- 
dimensional (ID) condensate of ^Li, loaded in a strongly elongated (nearly one-dimensional) "cigar-shaped" 
trap 0]. The use of the Feshbach resonance (FR) made it possible to keep the scattering length in the 
condensate negative, with a very small absolute value, ~ 0.1 nm. The weak self-attractive nonlinearity 
controlled by means of this technique was sufficient to create ID solitons, being a way below the collapse 
threshold in the cigar-shape configuration, thus securing the stability of the solitons. 

More generic for BEC is a positive scattering length, corresponding to repulsive interactions between atoms. 
In this case, bright solitons may be created as a result of the interplay of the intrinsic repulsion and periodic 
potential induced by an optical lattice (OL, i.e., the interference pattern created by counterpropagating 
beams illuminating the condensate). It was predicted [Hl^ that gap solitons (GSs) may emerge in bandgaps 
of the system's spectrum, since the interplay of a negative effective mass, appearing in a part of the gap, 
with the repulsive interaction is exactly what is needed to create a bright soliton. Theoretical models for 
GSs in BEC were reviewed in Ref. 3, and a rigorous stability analysis for them was developed in Ref. 0. 
Experimentally, a GS in the *^Rb condensate with a positive scattering length, loaded into a cigar-shaped 
trap equipped with a longitudinal OL, was for the first time created in Ref. 6] (the soliton was composed 
of- 1000 atoms). 

Binary mixtures of BECs are also available to experimental studies. Most typically, they contain two 
different hyperfine states of the same atomic species, such as ^^Rb and ^'^Na Q- BEC was also created in 
a heteronuclear mixture of '*^K and ^^Rb As mentioned above, the magnitude and sign of the scattering 
lengths of collisions between atoms in the same species may be altered, via the FR technique, by an external 
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spatially uniform magnetic 10] or optical [Tl| field. The scattering length of collisions between atoms 
belonging to different species may also be controlled by the magnetic field [ij ■ 

The latter possibility opens a way to create a binary mixture with the (natural) intra-species repulsion, 
while the sign of the inter-species interaction is switched to attraction. In recent preprints it was 
proposed to use this possibility to create symbiotic bright solitons: while each self-repulsive species cannot 
support a soliton by itself, the inter-species attraction opens a way to make two-component solitons. A 
somewhat similar possibility was earlier proposed in terms of a Bose-Fermi mixture, where the interaction 
between bosons is repulsive, but the bosons and fermions attract each other . Related to the latter setting, 
is a proposal to use attraction between fermions and bosons to built bosonic quantum dots for fermions (in 
particular, gap solitons in the BEC trapped in an OL may play the role of such quantum dots) [T^ . 

In the present work, the aim to construct 2D and ID solutions for two-component GSs in the most 
natural setting, when the inter-species interaction is repulsive. We will chiefly focus on the basic case, when 
intra-species interactions may be completely neglected, while the interplay of the OL potential and repulsion 
between the two species help to build GSs. The actual possibility to nullify the intra-species scattering length 
by means of the FR depends on the atomic species: as is known, it can be done in ^Li (see, e.g., Ref. [l Tj ) , 
while in ®^Rb loss effects grow close to the FR point. Another possibility is offered by spinor condensates, 
where the scattering lengths which determine collisions between atoms with the same or opposite values of 
the hyperfine spin, nip = ±1, can be represented, respectively, as |3| a = ao ± 02, the coefficients oq and 
ao accounting for the mean-field (spin-independent) and spin-exchange interactions between the atoms. In 
this case, the self-scattering length vanishes in the case of 02 = — oo- 

Besides that, the model with zero interaction inside each species and repulsion between them may also 
apply to a mixture of two ultra-cold Fermi gases |l5l Il6| . In this connection, it is relevant to mention that 
the scattering length of collisions between fermionic atoms (such as ^Li) may also be controlled by means 
of the FR jTvUTsj l. An effect of the intra-species repulsion will be briefly considered too, with a conclusion 
that it additionally stabilizes two-component GSs. 

The periodic OL potential gives rise to many bandgaps in the system's spectrum. In this work we con- 
centrate on the most fundamental situations, with the two components of the soliton belonging to two 
lowest-order gaps. This way, we will demonstrate intra-gap solitons, with both components sitting in either 
the first or second gap, and inter-gap solitons, with the components belonging to the different (first and 
second) gaps. 

The paper is structured as follows. The model is set in Section 2. An analytical variational approximation 
for 2D two-component solitons is presented in Section 3. Direct numerical results, that identify existence 
and stability regions of the intra- and inter-gap solitons, are reported in Section 3. The stability of the 
2D solitons is investigated in direct simulations, which reveal not only stable stationary solitons, but stable 
breathers too. Basic results for the ID version of the same two-component model are collected in Section 4; in 
particular, the stability of the ID soliton is identified via computation of eigenvalues for small perturbations 
(which is more difficult in the 2D case). The paper is concluded by Section 6, where we also give estimates 
for actual numbers of atoms in the solitons predicted in this work. 



In the mean-field approximation, the binary BEC at zero temperature is described by a system of two 
coupled Gross-Pitaevskii equations (GPEs) for the wave functions ^(X, Y,Z;T) and ^{X,Y, Z;T) of the 
two species p^: 



II. THE MODEL 
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where m is mass of of both species of atoms, Vq and n/k are the ampHtude and period of the OL potential, 
U{Z) is a potential accounting for the tight magnetic or optical confinement in the transverse direction, that 
makes the condensate effectively two-dimensional (squeezing it into a "pancake" shape |l9l|). while a and pa 
(with p > 0) are the scattering lengths of the inter-species and intra-species collisions. The equations do 
not include an external trapping potential in the {X, Y) plane, as we are interested in localized 2D states 
supported intrinsically by the interplay of the inter-species repulsion and OL potential. 

Assuming that the transverse trap gives rise to a ground-state wave function Xoi^) with the respective 
energy Eq, reduction of Eqs. J^l to a normalized system of effective two-dimensional equations follows 
the usual procedure, based on averaging in the Z direction and rescaling p^ . To this end, we define 
{kX,kY} = {x,y} , (^hk^/2rnjT = t, 2mVo/ (hk) = e {e measures the ratio of the height of the OL's 
potential barrier to the atom's recoil energy) , and 



The eventual form of the 2D equations is 



f'^°^XoiZ)dZ k 

Wx,y;t),cj^ix,y;t)}. (2) 



\ j^::xt>iz)dz 
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ifpt + ipxx + ipyy+e [cos(2x) + cos(2?/)] ip - (plV"!^ + |<^|^) V' = 0, 

(3) 

iipt + (pxx + (pyy+e [cos(2a;) + cos(2y)] (p ~ {p\p\'^ + IV'P) p = 0- 
In addition to e and p, control parameters of the normalized system are norms of both components, 

Ni = [ [ \'4j{x,y f dxdy, N2 ^ f f \ip{x,y)f dxdy, (4) 



which determine the respective numbers of atoms as per Eq. 
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Norms Q, together with the Hamiltonian of the system, are its dynamical invariants. As said above, we 
will be mostly dealing with the case of p = (no intra-species interactions), hence essential parameters are 
e and Ni^2 (generally, we will consider the asymmetric case, with iVi ^ N2)- 
Stationary solutions of Eqs. Q are looked for in the usual form, 

y, t) = u{x, y) exp(-i^ii), p{x, y, t) = v{x, y) exp{-ip,2t), (6) 

where /ii and p2 are real chemical potentials, and the real functions u(x,y) and v(x,y) are solutions of the 
equations 

Piu + u.j;x + Uyy + £ [cos(22:) + cos(2?;)] u — v^u = 0, 

(7) 

II2V + Vxx + Vyy + £ [cos(2a::) -|- cos(2?;)] v — u^v = 0. 

Linearization decouples Eqs. Q, hence each chemical potential must belong to one of bandgaps of the 2D 
linear operator, 



92 Q2 

^ ^ g;^ + q:^ + ^ [cos(2x) + cos(2y)] (8) 



(the spectrum of the operator is known, see, e.g., Ref. Q and Fig. below). Note that /ii and 112 may be 
placed in different gaps, which gives rise, as said above, to two-component inter-gap solitons, as opposite to 
ones of the intra-gap type, with jii and p2 belonging to the same bandgap. Below, we will demonstrate that 
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the model supports stable solitons of both types (the consideration will actually be limited to two lowest 
gaps, the intra-gap solitons being possible in both of them). 

It is relevant to mention that, in nonlinear optics, one-component ID spatial solitons sitting in higher 
gaps were experimentally observed in arrays of waveguides with the Kerr (cubic) nonlinearity |2(]| | . Then, 
2D spatial single-component lattice solitons with embedded vorticity, belonging to the second gap, were 
created in a photorefractive material equipped with a 2D photonic lattice ||21|]. Moreover, in the latter case, 
two-component solitons composed of a fundamental soliton in the first gap and a vortex soliton in the second 
gap were also created. The most essential difference of the photorefractive media from BEC is the saturable 
character of the photorefractive nonlinearity. 



III. VARIATIONAL APPROXIMATION FOR TWO-DIMENSIONAL SOLITONS 

It is known that ID and 2D solitons generated by the GPE with the lattice potential (in the 2D case, the 
lattice may be either two-dimensional or quasi-one-dimensional 23]) are naturally classified, in both cases of 
the attractive and repulsive (2^ interaction, as tightly-bound (TB, alias single-cell) and loosely-bound 

(LB, alias multi-cell) ones. The solitons of the TB and LB types are essentially localized in one or several 
cells of the OL potential, respectively. 

TB solitons, in the ID and 2D geometry alike (including the case of the 2D equation with the quasi-lD 
lattice), may be accurately predicted by the variational approximation (VA) E| (this approximation 
was first applied to the GPE in Refs. [2g; a general review of the technique can be found in Ref. 27]). For 
LB solitons, the VA may provide an adequate approximation for the soliton's central lobe only, but not for 
the slowly decaying oscillatory tails 28]. 

We start the analysis with application of the VA to Eqs. {Tj). The Lagrangian from which the equations 
are derived is L = / J Cdxdy, with the density 

C = niu^ + /i2i;^ -ul-ul-vl-vl 

+e [cos(2a;) -f cos(2y)] {u^ + v^) - u'^v'^. (9) 

We approximate the TB solitons by a simple isotropic Gaussian ansatz, 

u(x,y) = Aexp(^-^^±^^ , v{x,y)^Bcxp(^-^^y (10) 
The substitution of the ansatz in Eq. ^ and integration yield an effective Lagrangian, 



N^N■? 

■ A^iA^i+M2^2, (11) 



TT (a2 + 52) 



where the norms of the components, Ni — ttA^o^ and N2 = irB^b^, are obtained by the substitution of the 
ansatz (|1U|) in Eqs. Q. Below, we use, instead of A^i and N2, the total and relative norms, 

N = Ni + N2, Nr = N1/N2 (12) 

(we define Ni and N2 as the smaller and larger norms, respectively, hence Nr < 1). 

The variational equations, dL/dNi 2 = dL/da = dL/db — 0, applied to Lagrangian Hll|l . yield the 
following relations that determine parameters of the soliton within the framework of the VA: 

a'^ I - 2eb^e-''' 

= 7rfl + — j {l + Nr)(2sa-^e-''' -ij, (14) 
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l)ee-''\ (15) 
l)ee-^\ (16) 

Equations \VA\ and H14(l give rise to necessary conditions for the existence of the sohton, 1 — 2ea'*e^" < 
and 1 — 2elAe^^ < 0, from which, in turn, it follows that that the GS does not exist unless the OL strength 
e exceeds a minimum (threshold) value, Emin > « 0.92. The same condition for the existence of GSs 
was predicted by the VA in the single-component 2D model [28|. 

The VA does not include information about the location of bandgaps in the system's spectrum, therefore 
it is necessary to check whether the GSs predicted by the VA fall into the bandgaps. This is shown in 
Fig. ^ which demonstrates that the variational equations H13|) . H14() and H15() . (|16|l for the symmetric GSs 
{Hi = H2 = fJ,) predict the boundary of the first finite gap (the lower bold dashed line) with surprisingly good 
accuracy. The accuracy may be explained by the fact that the Gaussian ansatz (llUf) provides for a good fit 
to the actual shape of the GS in the first finite gap (this, in particular, means that the VA correctly predicts 
the absence of solitons in the semi-infinite gap - the lowest one in Fig. ^ which extends to /i ^ — oo). On 
the other hand, the correlation between the other VA-predicted soliton-existence boundary (the upper bold 
dashed line in Fig. ^ and the exact border of the ba ndg aps is very crude, due to the fact that the GSs in 
higher gaps are very different from the simple ansatz j28l | . 

For the description of experimentally relevant situations, Fig. [5]displays the predicted GS existence regions 
in the plane of the total and relative soliton norms, {N, Nr) [see the definitions in Eqs. (|12l) ]. which includes 
the general asymmetric solutions with A^^ < 1- False parts of the existence regions, that have either /ii or 
/i2 falling into a Bloch band (rather than into a gap), are excluded in panels (b)-(d). 

We note that, for relatively small e (however, it must exceed the above-mentioned threshold value Emin 
necessary for the existence of the solitons), the VA predicts only intra-gap solutions, with both fii^2 lying in 
the first finite gap. This prediction is correct, as for these values of e the spectrum of the two-dimensional 
GPE with the OL potential has only one finite gap, see Fig. ^ For larger e, the VA predicts intra-gap 
solutions in the second gap, as well as inter-gap solitons, with the components belonging to the different 
finite bandgaps, first and second. With the further increase of e, the regions of the intra-gap solutions in the 
first and second bandgaps and inter-gap solitons grow and overlap. 



Ml = j + 2(a^ + b^- 



IV. NUMERICAL RESULTS FOR TWO-DIMENSIONAL SOLITONS 



A. Families of intra- and inter-gap solitons 

Comparison with numerically found stationary solutions for the GSs demonstrates that (as mentioned 
above) the applicability of the VA is indeed limited to strongly confined TB solitons. In most cases, the 
solitons with both components belonging to the first finite gap appertain to this type, and they are predicted 
by the VA very accurately, as shown in Fig. |21 Note, however, that strongly asymmetric intra-gap solitons 
cannot be found numerically in the first bandgap for values of the relative norm smaller than (Nr)^-^^^^ ~ 0.05. 
Nevertheless, the VA predicts asymmetric solitons for Nr < {Nr)^^^, up to Nr = 0. The latter is an obvious 
artifact of the VA, as in the case of Nr = one component is empty (^p = 0), which makes the remaining 
equation linear, hence it cannot give rise to any soliton. 

On the contrary to the situation for the intra-gap solitons belonging to the first finite bandgap, the 
prediction produced by the VA for the solitons with both components sitting in the second bandgap are 
completely wrong: as seen in the top panel of Fig. 13 the formally predicted family of the intra-gap solitons 
in the second gap has no numerically found counterpart. 

Increase of the total norm N pushes the solution to the higher (second) bandgap. Figure 01 compares the 
numerical results and VA for such a situation (with five times as large as in Fig. We observe that, 
while the VA is good for the symmetric solitons {Nr = 1), it fails for more general asymmetric solutions - 
both intra-gap solitons belonging to the second gap, and inter-gap ones. 

Figure 01 also shows that, as the asymmetry coefficient Nr gets smaller, the soliton's component with the 
smaller norm loses its single-peaked shape, see panels (c)-(e). We stress that no intra-gap soliton belonging 
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to the second gap could be found (numerically) for Nr < 0.75. The latter is seen in Fig. ^ where the 
branches corresponding to the intra-gap solution family terminate at ~ 0.75. 

The findings presented in Fig. ^include inter-gap solitons [in panels (d) and (e)], with ^2 in the first gap 
and /^i in the second. As seen from panel (a) in the figure, these solitons could not be predicted by the 
VA. This is explained by the fact that at least one component of the inter-gap soliton has a loosely-bound 
shape, hence the Gaussian ansatz (|10|l is irrelevant for it. A typical example of the inter-gap soliton built 
as a TB-LB bound pair is shown in Fig. [S] for a still larger total norm, N — 70. As well as in the above 
examples, the TB and LB components reside in the lower and upper gaps, respectively. 

Intra-gap solitons may also be TB-LB bound states. Figure demonstrates how a symmetric TB soliton 
(with Nr = 1) transforms into a TB-LB pair with the decrease of N^, this time (on the contrary to the above 
examples of the inter-gap TB-LB bound states) the LB component having a larger norm, which is, in fact, 
a rule for asymmetric intra-gap solitons, as evidenced by our numerical results obtained with other values of 
the parameters. 

With even larger norms (for instance, N — 150), intra-gap solitons in the form of LB-LB bound pairs with 
very slowly decaying tails were found too, for Nr > 0.85, i.e., they are nearly symmetric states. However, 
such solitons are unstable. 



B. Global existence and stability diagrams for the solitons 

Results of systematic investigation of the existence and stability of the two-component two-dimensional 
GSs arc summarized in Fig. |7| which is a typical example for the weak OL, with e — 2, that supports only 
one finite bandgap in the spectrum of the 2D model, and in Fig. |H| which represents relatively strong OLs 
(it has £ — 10, which admits two distinct finite bandgaps). 

Figure [3 shows that the entire gap is populated with solitons which are stable, except when both chemical 
potentials /zi and 112 are close to either the lower or upper bandgap edge. The stability was verified by direct 
simulations of the underlying GPEs ||2J), with an initial perturbation imposed on the soliton by dislocating 
centers of its two components (simulations were performed by means of the split-step algorithm combined 
with the 2D fast Fourier transform). The dislocation leads to oscillatory dynamics, and GSs were classified as 
stable ones if they would oscillate near the initial shape. A caveat is that such a definition of the stability does 
not make a clear distinction between stable stationary solitons and stable breathers with a small amplitude 
of internal vibrations; however, the objects of both types are, as a matter of fact, experimentally equivalent 
localized states in the BEG. 

Figure |S| shows that the stability pattern is more complex for a stronger OL, with e — 10. In particular, 
the inter-gap soliton, which is possible in this case, may be stable when the chemical potential /ii of the 
component with a smaller norm, which appertains to the first (lower) bandgap, is sufficiently close to the 
lower edge of the gap. It is noteworthy that the stability region of the inter-gap solitons is found at values 
of the total norm N at which intra-gap solitons, with both components sitting in the first bandgap, do not 
exist, i.e., there is no overlap between these two types of the stable GSs . 

Further simulations demonstrate that unstable inter-gap solitons are either completely destroyed or evolve 
into symmetric solitons with both components belonging to the first gap. An example of such evolution is 
displayed in Fig. ^ In this case, the LB component of the TB-LB pair emits radiation until only the central 
lobe remains in it, and the whole structure turns into a symmetric intra-gap soliton of the TB type. 

A specific instability mode was observed in unstable inter-gap solitons, as well as in intra-gap ones with 
both components belonging to the second gap, in the case when the unperturbed soliton features a single 
peak in one component and a double (split-tip) peak in the other. In this configuration, the instability 
triggers oscillatory dynamics with the single peak jumping irregularly between positions close to the two 
side peaks of the mate component. An example of this instability is displayed in Fig. 1101 Further evolution 
leads to complete destruction of both components in this case. 

Finally, simulations of the full model © with p 0, which includes the self-repulsion in each component, 
reveal a clear trend to stabilization of the two-component GSs as p increases (of course, the model with 
p > gives rise to ordinary single-component GSs too). For example, the unstable symmetric soliton with 
pi — p,2 = 1.2, s = 2, and p = (recall the system has a single finite bandgap in this case) becomes stable 
as the self- repulsion coefficient increases to Pmin ~ 2/3. This stabilization mechanism will be considered in 
more detail below in the ID version of the model. 
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V. ONE-DIMENSIONAL SOLITONS 



The ID variant of the model corresponds to a "cigar-shaped" binary condensate, tightly confined in the 
transverse plane, with a one-dimensional OL created in the longitudinal direction, x. The accordingly 
modified ID version of the normalized equations (PJ is 

i^l^t + ^xx+£C0s{2x)il:~{p\il)\'^ + \lp\^)iI) = 0, 

(17) 

iipt + ipxx+ £C0s{2x)ip ~ [p\ip\'^ + \'tj)\'^) ip = 0. 
The conserved norms of the scaled wave functions are 

/ + 00 
[\i^{x)\\\^{x)\'\dx, (18) 

which are related to the numbers of atoms as fohows: 

cf. Eq. ©. Here vr/fc is, as above, the OL period, and xo(^) is the ground-state wave function of the tight 
confining potential, R being the radial coordinate in the transverse plane. 

A principal difference of the ID model from its 2D counterpart is that, at any finite e, the ID version 
of the operator JSJ , which is the same as in the Mathieu equation, gives rise to an infinite system of finite 
bandgaps (recall that the 2D operator generates a single finite gap for small e, two gaps for larger e, etc.). 

Stationary soliton solutions to Eqs. (|17|l are again looked for in the form of Eqs. © (with the coordinate 
y dropped). In contrast to the above analysis of the 2D model, in the present case we determine the stability 
of solitons in a rigorous way, from linearized equations for small perturbations about the stationary soliton 
(however, in all the cases when GSs were predicted to be stable in this sense, their stability was also verified 
in direct simulations). The application of the rigorous approach to the 2D model is to be presented elsewhere, 
as it is a technically involved problem. 

The perturbed solutions are taken as 



tp = [u{x) + ^pi{x)e 

(20) 

= [vix) + ipi{x)er'^'] e^'f^^^ 

where ipi and ipi are eigenmodes of infinitesimal perturbations and A is the respective eigenfrequency, the 
instability corresponding to having Im A > 0. The substitution of this in Eqs. (|17|l and linearization lead to 
the equations 

^ (f \ 

P2 + j ^1 + e cos{2x)ijji - [v'^{x)^i + u{x)v{x) (^i + ^i)] 

-pu^ {x)i2i;i+ri) = AV-i, 

^A^i + '^i + e cos{2x)(t)i - [u'^{x)<j)i + u{x)v{x) {tpi + -0^)] 

-pv\x){2<j)i + cl)l) = A01, 

which can be solved by means of known numerical methods, to yield a full spectrum of the eigenfrequencies A 
(we used the Matlab eigenvalue-finding routine; it is based on approximating the ordinary differential equa- 
tions by a system of linear homogeneous algebraic equations, computing the determinant of the corresponding 
matrix, and equating it to zero, which eventually leads to an equation for A). 

First of all, we present results for the most fundamental case of p = 0, when only two-component GSs are 
possible, as well as in the 2D model considered above. Fixing the OL strength e, in Figs. Illf a') and (b) we 
display the total and relative norms, N and Nr for the family of ID two-component GSs, as functions of 
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the two chemical potentials, fii and /12, for a case when both belong to the first finite bandgap (i.e., for the 
family of intra-gap solitons) [the norms N and Nr are defined as in Eqs. H12() . with A^i and N2 taken as per 
Eq. H18|) ]. The stability of the same GS family is presented in panel (c) of Fig. ^2 which shows the largest 
value of Im A, i.e., the instability growth rate, as a function of ni and fi2 [the border between stable and 
unstable solitons is indicated in Figs. Illf a) and (b) too]. If the GS of this type is unstable, its instability 
is oscillatory (i.e., Im A > comes along with Re A 7^ 0). Typically, the instability does not destroy the 
soliton, but rather transforms it into quite a stable breather, see an example in Fig. IllT d). 

Inter-gap solitons, built of components belonging to the fist and second finite bandgaps, have also been 
investigated in detail in the ID model. Characteristics of this solution family are presented in Fig. E| As 
seen from panel (c) of the figure, the stability area of the inter-gap solitons is essentially smaller than in the 
case of their intra-gap counterparts, cf. Fig. Hir e). Note that panels (a), (b) and (c) in Fig. E|do not show 
a Bloch band that separates the two gaps (unlike Fig. |Slin the 2D model), as the ranges of values on the 
axes of ^1 and ^2 in these panels cover, respectively, only the first and second gaps. 

It is noteworthy that, as well as in the 2D model, the inter-gap solitons are bound states of TB and LB 
components belonging to the first and second gaps, respectively. An example of a stable inter-gap soliton 
that clearly demonstrates this structure is shown in Fig. I12r d) [cf. Figs. Eld) and (e) and Fig. |S1 which 
display cross sections of intra-gap solitons in the 2D model]. 

If an inter-gap soliton is unstable, its instability again has the oscillatory character. However, the action of 
the instability on the inter-gap soliton is more destructive than it was in the case of its intra-gap counterpart: 
instead of transforming the stationary soliton into a well- localized breather [see Fig. Illf d)]. the instability 
triggers much more violent evolution, as illustrated by a typical example in Fig. 1131 

One-dimensional intra-gap solitons belonging to the second gap were constructed and investigated too, 
with a conclusion that they all are unstable, although the instability growth rate may sometimes be very 
small. An example illustrating the instability of this species of the GS is displayed in Fig. E| This conclusion 
seems to be in contrast with results reported above for the 2D model, where a small stability area for the 
intra-gap solitons appertaining to the second finite bandgap was found, see Fig. |S1 However, the stability of 
the 2D solitons was identified not via eigenfrequencies of small perturbations, but rather by means of direct 
simulations, and the numerical stability test always turned the stationary soliton into a weakly excited state 
(a breather with a small amplitude of intrinsic vibrations). Therefore (as it was said above), in the 2D case 
we, strictly speaking, cannot tell a difference between a weakly unstable stationary soliton and a weakly 
excited stable breather into which the dynamical evolution may transform the unstable soliton. Thus, it 
may happen that what was identified as stable 2D intra-gap solitons sitting in the second finite bandgap are, 
actually, stable breathers with a small vibration amplitude, similar to what is observed in Fig. Illf d'). 

We have also investigated in some detail the effect of the self-repulsion term in the full ID system IjlTI) . As 
was already mentioned above in connection to the 2D model, the increase of the self-repulsion coefficient p 
leads to stabilization of the GSs. A similar effect in the ID setting is clearly seen in Fig. ^1 It demonstrates 
that a stable symmetric intra-gap soliton belonging to the second finite bandgap, which, as said above, is 
always unstable in the ID model with p = 0, becomes stable if p exceeds a minimum (threshold) value, 
which, in this case, is pmin ~ 0.9. 

The above stability analysis was performed within the framework of the GPE, i.e., in the mean-field 
approximation at zero temperature. A physically important issue is stability of the same solitons against 
quantum and thermal fluctuations beyond the framework of the mean-field theory. The interest to this issue 
is enhanced by experimental observation of quite strong effective re-thermalization the ultra-cold condensate 
under the action of the OL potential |29j] . To extend the analysis in this direction, one may use a system of 
self-consistent time-dependent Hartree-Fock-Bogoliubov (TDHFB) equations, built around the corresponding 
Gross-Pitaevskii equation(s) (30| . This approach was recently used to demonstrate possible instability of an 
ordinary ID soliton (not of the gap ty pe) in the BEG with attractive interactions, which is completely stable 
in the framework of the GPE proper |3l|. In that case, the instability splits the soliton into two segments. 

A properly modified system of the TDHFB equations (including the OL potential and repulsive inter- 
species interaction) can also be used for the investigation of the extended stability of the two-component 
GSs studied in this work. We have performed a preliminary analysis along these lines, and concluded that 
those ID solitons which are stable within the framework of Eqs. (II 7|) are also stable (at least, in most cases) 
against fluctuations governed by the TDHFB equations. A systematic consideration of this issue, especially 
for the 2D model, requires a considerable amount of additional work and will be reported elsewhere. 
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VI. CONCLUSION 



In this work, we have introduced a model of a binary BEC with intrinsic inter- and intra-species repulsion 
(positive scattering lengths), which is loaded into the periodic optical-lattice potential. Both two- and one- 
dimensional versions of the model were considered. We focused on the most fundamental case (different from 
previously studied models) , with repulsion between the two species and zero intra-species interaction, which 
can be achieved by means of the Feshbach-resonance (FR) technique, or in a spinor condensate. The same 
system may also model a mixture of two mutually repulsive fermionic species. 

The main problem was the existence and stability of gap solitons (GSs) supported by the interplay of 
the inter-species repulsion and periodic potential. Two-component GSs were looked for by means of the 
variational approximation (VA) and in the numerical form. It was found that the VA provides for good 
accuracy in the case when the 2D soliton is a bound state of two tightly-bound components, each being 
essentially confined around one cell of the periodic potential. Such a GS structure dominates in the case 
when both components belong to the first finite bandgap of the system's spectrum (the intra-gap soliton). 
In fact, only this type of the GS is possible in a weak 2D lattice potential. On the contrary, inter-gap 
solitons, and intra-gap ones residing in the second finite bandgap (both types are possible in a stronger 
lattice potential, with the barrier height essentially larger than the recoil energy) are, typically, bound states 
of tightly- and loosely-bound components. For such structures, the VA is irrelevant, but general results were 
obtained in the numerical form, which made it possible to identify the existence and stability regions for the 
inter- and intra-gap solitons in both the 2D and ID models. In the 2D case, the stability was tested in direct 
simulations, while in the ID model the stability was identified in a rigorous way, through the computation of 
eigenfrequencies of small perturbations (the results were also verified by direct simulations). In the case when 
the ID intra-gap soliton belonging to the first gap is weakly unstable, it evolves into a stable breather with 
a small amplitude of intrinsic vibrations. In contrast to this, if the 2D solitons in the first gap are unstable, 
they are completely destroyed by the instability. The same pertains to unstable ID and 2D solitons of other 
types, such as inter-gap solitons, and intra-gap ones belonging to the second finite bandgap. It was also 
shown that introduction of the intra-species repulsion, in addition to the repulsion between the components, 
leads to further stabilization of solitons. In particular, some originally unstable types (such as ID intra-gap 
solitons in the second bandgap) may be made stable, provided that the self-repulsion coefficient exceeds a 
certain minimum value. Preliminary analysis shows that the two-component GSs introduced in this work 
are stable too against fluctuations obeying the Hartree-Fock-Bogoliubov equations. 

The actual number of atoms in the 2D solitons considered above, which is their most important physical 
characteristic, can be easily estimated. Undoing renormalizations which led to the 2D Gross-Pitaevskii 
equations in the form of Eqs. (O and making use of Eq. jSJ, one can derive the following relations between 
the density of atoms in physical units, nphys, and the rescaled one IV'Pj and between the number of atoms 
and the soliton's norm: 

where A is the period of the optical lattice, the scattering length of the inter-species collisions, and Iz the 
size of the condensate in the transverse (z) direction (note that the relation between iVphys and N does not 
contain A). Assuming experimentally reasonable values, ~ 1 nm and 1^ — 2 /.tm (the latter pertains to 
a "p ancake-shaped" condensate trapped between two difi^raction-limited blue-detuned repelling laser sheets 
[l9|'). we conclude that typical values, TV ~ 20 and N ~ 100, for the stable intra-gap and inter-gap solitons, 
respectively [see Figs. EJa) and|HIa)], lead to the following estimates for the respective numbers of atoms: 
Mntra ~ 5, 000, TVintor ~ 25, 000. If ncccssary, this number may be made at least an order of magnitude 
larger, reducing by means of the FR technique IT^. It is relevant to mention that, in the first observation 
of a GS in BEC ^6], the number of atoms ~ 1,000 in the soliton's central lobe was quite sufficient for the 
experiment. 

For the effectively ID condensate in a cigar-shaped trap Eq. H19|) leads to the following equation 
replacing the second relation in Eq. (|^ . 

N,.y. = ^N, (22) 



where A is an effective transverse area of the trap [note that this relation, unlike its 2D counterpart in Eqs. 
H21|l . explicitly contains the optical-lattice period A]. Taking the same as above, ^ 1 nm, a typical A ~ 1 
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/^m, and a reasonable value of A ~ 30 /im^ (it corresponds to the effective trap's radius ~ 3 /ini), we conclude 
that, with the characteristic value of the norm for the stable intra- and inter-gap ID solitons, iV ~ 10 in the 
normalized units [see Figs. Illf a) and ll2r a)]. Eq. 1)22(1 yields an estimate ~ 10^ for the actual number of 
atoms in the effectively ID soliton. 
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FIG. 1: The region of the existence of symmetric gap sohtons (the area between the two bold dashed 
lines), as predicted by the variational approximation based on the Gaussian ansatz H1U|) . and the ex- 
act bandgap structure in the two-dimensional system. Shaded and unshaded regions are, respectively, 
the Bloch bands (where solitons cannot exist) and gaps (where solitons are possible), the numbers 
1, 2, 3 being numbers of the finite bandgaps. The lower unshaded region is the semi-infinite gap. 



13 













450 




\ 

\ 






400 




\ 






350 




\ 

N 






300 




e = 10 






Z 250 










200 


-\ 

\ 








150 




8 = 4 






100 










50 




8 = 2 








( 












0.2 0.4 


0.6 


0.8 



(a) 



8 = 2 



Gap 1 



0.4 0.6 0.8 



(b) 



8 = 4 




Gap 1 



0.4 0.6 

N 



(c) 




FIG. 2: (a) The existence regions for the gap sohtons in the {N,Nr) plane, as predicted (below the 
respective dashed Unes) by the variational approximation for different values of the strength e of the 
optical lattice. Panels (b)-(d), which pertain, respectively, to e = 2, 4, and 10, show that shaded 
portions of the predicted existence regions must be excluded, as in these areas cither chemical poten- 
tial, /Lti or fj,2, falls into a Bloch band. The remaining white portions of the existence region in each 
panel are true ones, with both /^i and ^2 located in one of the two lowest bandgaps. In panels (b) 
and (c), the bandgap(s) to which the chemical potentials belong are indicated. In panel (d), details 
for different gaps are not shown, because of a complex picture of areas corresponding to each gap. 
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FIG. 3: (Color online) Comparison of the variational approximation and numerical solutions for two- 
dimensional gap solitons in the case of e = 10, with the total norm N = 10. The soliton fam- 
ilies, as predicted by the VA and found in the numerical form, are shown by dashed and solid 
lines, respectively, in panel (a). In this panel, areas occupied by the Bloch bands (where soli- 
tons cannot exist), are shaded. Small circles (o) and crosses (x) designate, respectively, variational 
and numerical solutions that are displayed, as examples, in the remaining part of the figure: (b) 
Nr = 1 [a. symmetric soliton); (c) = 0.5; (d) = 0.05. In these panels, and in exam- 
ples of the two-dimensional solitons presented below, their cross-sections along y — Q are displayed. 
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FIG. 4: (Color online) The same as in the previous figure, but for N = 50. Panels (b) and (c): ex- 
amples of the intra-gap solitons, for A'^^ = 1 and Nr = 0.8, respectively, which belong to the second 
bandgap. Panels (d) and (e): examples of inter-gap solitons, with A'^^ = 1 and Nr = 0.1, respec- 
tively [in both these cases, the loosely-bound component belongs to the second (upper) bandgap, and 
in the latter case {Nr = 0.1) the larger norm is in the first (lower) bandgap]. Note that the inter- 
gap soliton with Nr = 1 is not a symmetric one, even if the norms of its two components are equal. 
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FIG. 5: An example of an inter-gap soliton, built as a bound state of tightly- 
and loosely bound components with equal norms (Nr = 1), for e = 4 and 
N = 70. The corresponding chemical potentials are = —1.64 and /X2 = 3.3. 
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FIG. 6: (Color online) Stable intra-gap solitons numerically found for e = 2 (in this case, the sys- 
tem's spectrum contains only the first finite bandgap), with the total norm N = 20. (a) The 
family of the gap-soliton solutions. (b) A tightly-bound symmetric soliton with Nr — 1. (c) 
An example of a strongly asymmetric soliton, with N,. — 0.1, in the form a bound state of 
loosely- and tightly bound components, with the larger norm sitting in the loosely-bound component. 
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FIG. 7: (Color online) The global existence and stability diagram for the two-component intra-gap soli- 
tons in the weak two-dimensional lattice potential, with £ = 2 (the spectrum of the 2D Gross-Pitaevskii 
equation contains only one finite bandgap in this case). The total norm of the soliton, N = Ni + N2, 
(a), and the relative norm, Nr = N1/N2, (b), are shown vs. the chemical potentials of the two compo- 
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FIG. 9: (Color online) Evolution of an unstable two-dimensional inter-gap soliton with /ii = 
—3.7, fi2 = —5.6 is shown in the cross-section along the y — axis for e = 
10. The loosely-bound component reduces its norm through emission of radiation. Even- 
tually, the solution evolves into a symmetric intra-gap soliton of the tightly-bound type. 
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FIG. 10: (Color online) The onset of the "jumping" instability of the two-dimensional intra-gap soli- 
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FIG. 11: (Color online) A typical example of a family of two-component intra-gap solitons in the one- 
dimensional model, with both components belonging to the first finite bandgap. In this case, e = 8 (and 
p = 0). Panels (a), (b), and (c) show, respectively, the total norm N = Ni + N2, relative norm Nr — N1/N2, 
and instability growth rate, Im A. Panel (d) displays an example of the transformation of an unstable 
symmetric intra-gap soliton (with fJ-i = fJ-2 ~ 0) into a stable breather, due to the oscillatory instability. 
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FIG. 12: (Color online) Panels (a), (b), and (c) show the same characteristics as in Fig. ^2 for a typical 
family of two-component inter-gap solitons in the one-dimensional model, with one component belonging 
to the first finite bandgap, and the other - to the second. In this case, e = 8 and p — 0. Panel (d) 
displays a generic example of a stable inter-gap soliton, with /ii = —3, A'^i = 11.7, and 112 = 3, -/V2 = 0.75. 
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FIG. 13: (Color online) An example of the evolution of an unstable inter-gap soliton in the one-dimensional 
model with e — 8, p = and /ii = —3, fi2 — 4. Panels (a) and (b) show the development of 
the instability in the components that originally belong to the first and second finite gaps, respectively. 




FIG. 14: (Color online) Instability of a symmetric intra-gap soliton appertaining to the second fi- 
nite bandgap in the one-dimensional model. In this case, e = 8, p — and p = 4. 
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FIG. 15: (Color online) The instability growth rate of a symmetric soliton in the second fi- 
nite bandgap vs. the self- repulsion coefficient p [see Eqs. (I17|l ]. Except for p, parameters are 
the same as in the previous figure. The soliton is stable in the region of p > Pmin ~ 0.8. 



